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Abstract 


The inverse problem of estimating the open channel flow roughness is solved using 
an embedded optimization model. Measured data for flow depths and discharges 
at different locations and times are used as inputs to the optimization model. The 
nonlinear optimization model embeds the finite-difference approximations of the 
governing equations for unsteady flow in Open-Channel as equality constraints. 
The Sequential Quadratic Programming Algorithm is used to solve the Optimiza- 
tion model. The performance of the proposed parameter estimation model is 
evaluated for different scenarios of data availability and noise in flow measurement 
data. Solution results for illustrative problems indicate the potential applicability 
of the proposed model. 



Contents 


Abstract iv 

Contents v 

List of Figures vii 

List of Tables viii 

Nomenclature ix 

1 Introduction 1 

1.1 Overview and Empirical Methods 1 

1.2 Literature Review 3 

1.3 Objective of the Study 6 

1.4 Organization of the Thesis 6 

2 Methodology 8 

2.1 Methodology 8 

2.2 Governing Equations 9 

2.3 Simulation Model 9 



CONTENTS 


vi 


2.4 Formulation, of the Inverse Problem 13 

2.4.1 Objective function 14 

2.4.2 Linear constraints 14 

2.4.3 Non linear constraints 14 

2.4.4 Optimization algorithm 16 

3 Results and Discussion 21 

3.1 Performance Evaluation 21 

3.1.1 Single channel problem 22 

3.1.2 Multiple reach problem 25 


4 Summary And Conclusion 


35 




List of Figures 


2.1 Flow Chart for Double Sweep Method of Solution 12 

3.1 Log Pearson Typelll Hydrograph 23 

3.2 Estimated and observed stage hydrographs at X= 20Km (Sce- 
nario 1) 24 

3.3 Estimated and observed stage hydrographs at the end of third 

reach (Scenario 9) 27 

3.4 Schematic layout of the multiple channel system 29 

3.5 Estimated and observed stage hydrographs 6000m downstream 

of junction point in channel 3 (Scenario 12) 29 

3.6 Effect of a on the parameter estimation 


34 



List of Tables 


3.1 Grid size effect on the parameter estimation 24 

3.2 Effect of initial estimate on the performance of the model .... 25 

3.3 Grid size effect: Multiple reach single channel problem 26 

3.4 Grid size effect ‘.Multiple channel problem 28 

3.5 Performance evaluation of the model for the multiple channel 

problem 31 

3.6 Estimated values of “n” for an a = 0.1 33 

3.7 Estimated values of “n” for a = 0.2 33 



Nomenclature 


ix 


N omenclatur e 


n 

h 

Q 

U 

s 0 

s f 

k 

9 

t 

x 

Ax 

At 

f 

3 

k 

9 

a 

N 

P 

x(o) 

x(m) 

A 

€* 


= Manning’s coefficient. 

= Depth of flow in the channel. 

= Dischage per unit width of the channel. 

= Velocity of flow in the channel. 

= Bottom slope of the channel. 

= Friction slope. 

= Conveyance of the channel. 

= Accleration due to gravity. 

= Time. 

= Longitudinal distance in metres. 

— - Grid spacing 
= Computational time step. 

= Any variable sj,q or h 
= Grid point in space. 

= Time level. 

= Weighting coefficient in Priessmann scheme. 

= Induced error percent. 

= Total number of nodes. 

= Total number of time steps. 

= Simulated value of depth or discharge from optimizer. 

= Observed value of depth or discharge. 

= Simulated noise free measurement data 
= Erroneous measurement data 

= A random error term sampled from a uniform distribution 



Chapter 1 


Introduction 


1.1 Overview and Empirical Methods 


An accurate estimation of Manning’s roughness coefficient is of primary im- 
portance in any study involving open channel flows such as: 

• Flood estimation, routing and damage mitigation, 

• Optimal design and operation of canal and irrigation systems, 

• Estimation of aggradation and degradation of rivers due to natural and 
man made causes, 

• Utilization of surface water resources and other river engineering studies. 

Estimation of Manning’s coefficient for a natural channel is not a trivial 
task as it depends on several factors like surface roughness, vegetation, channel 
irregularity, obstructions, channel alignment, sedimentation and scouring. Many 
empirical procedures (French, 1985) have been suggested in the past for esti- 
mating the value of “n” ,the Manning’s coefficient. The soil conservation service 
method (SCS)(Urquhart,1975) for estimating “n” involves selection of a basic 
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value of “n” corresponding to uniform, straight and regular channel in native 
material and then modifying its value by adding correction factors. The cor- 
rection factors are determined by a heuristic method by consideration of some 
of the above mentioned factors. 

A second method for estimating “n” utilizes of standard tables. Chow 
(1959) presented an extensive table containing a minimum, normal, and a 
maximum value of ”n” for various types of channels. This method gives a 
satisfactory estimate for man made channels, but is not suitable for natural 
channels. According to U.S Geological Survey, the photographs of channels 
of known resistance along with the knowledge of hydraulic and geometric pa- 
rameters which define the channel for a specified flow rate can be useful in 
estimating the resistance coefficient (Barnes, 1967). The U.S Geological sur- 
vey also maintains a program which trains the engineers in the estimation of 
“n” . However, this is not a formal procedure and the accuracy of estimation 
depends on the past experience of the user. 


From a theoretical viewpoint, value of the resistance coefficient can be 
estimated from velocity and depth measurements which depend upon rough- 
ness. Chow(1959), suggested that Manning’s “n”can be estimated with the 
help of velocity and depth measurements at a cross-section using an empirical 
formula given as follows. 

n = j£ziAA (u) 


6.78(z + 0.95) 


where 


X = U 0 . 2 /U 0 .s (1.2) 

f/o .2 = velocity at 2/10 th of depth 
U 0 8 = velocity at 8/10 i/l of depth. 
h = depth of flow. 

The above equations (1.1) and(1.2 ) are applicable only when the channel is 
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wide and velocity measurements are made for uniform flow conditions. 

This idea of estimating Manning’s “n” from field measurements of depth 
and discharge can be used for non-uniform flows conditions in irregular chan- 
nels. The Mannings “n” for the channel may be determined by a trial and error 
procedure such that the difference between simulated steady, nonuniform water 
surface profile using standard step method (Chow, 1959 and Subramanya,1984) 
and observed water surface profile is minimum. The above problem can be bet- 
ter formulated considering the following two aspects. Firstly, a formal search 
technique (Rao,1984) may be used which avoids trial and error, lends formality 
to calibration procedure, decreases arbitariness and gives likely error in the 
estimated value. Secondly, the best solution can be expected when abundance 
of data exist, which fortunately is the case with unsteady flow. In the present 
study a formal approach is developed for estimating the Manning’s ”n” from 
measured transient flow data. The developed methodology is based on non- 
linear optimization technique with the governing flow equations embedded as 
constraints. 


1.2 Literature Review 

The determination of parameters for a given system input and output and 
subject to boundary conditions is called an inverse problem. It is possible 
to estimate Manning’s “n” from field measurements of discharge and depth 
using Optimization techniques. Parameter estimation techniques using un- 
steady flow data are more popular because large number of data can be ob- 
tained from a limited number of gauging stations. Significant work has been 
reported in the area of ground water engineering relating to parameter esti- 
mation. Solution techniques for aquifer parameter identification include Quasi 
linearization used by Yeh and Tauxe (1971), Linear programming used by Klei- 
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necke(1971), Multiple objective Linear Programming used by Ne umann (1 973) 
and Marquaardt’s non-linear estimation algorithm used by Garay et al. (1976). 
However, the literature relating to identification of parameters in unsteady 
open channel flows is sparse. 

Becker and Yeh (1972) proposed a methodology for parameter estimation 
in unsteady open channel flow using influence coefficient approach. Friction 
slope and exponent of hydraulic radius in empirical friction slope relation were 
considered as the parameters to be identified. The methodology involves min- 
imizing the objective function which is the sum of the squares of error between 
observed data and simulated numerical solution of the equations based on 
some estimate of the parameters. An initial estimate of parameters is first 
specified and then these parameters are modified based on the evaluated ob- 
jective function and influence coefficient matrix at each iteration. Becker and 
Yeh (1973a) extended this model for identification of multiple reach channel 
parameters. The parameters estimated were bed slope and Manning’s ”n” at 
every reach. The numerical solution of governing equations in the above mod- 
els was based on explicit finite difference scheme method which requires a very 
small computational time step value. Also, the influence coefficient algorithm 
becomes tedious when used for estimating a large number of parameters. Yeh 
and Becker (1973) used Simplex linear programming algorithm with influence 
coefficient approach for estimating the unsteady open channel flow parame- 
ters, Manning’s ”n”, exponent of hydraulic radius and the lateral inflow. The 
minmax criterion (minimizing the maximum of absolute values of errors) was 
taken as the objective function. 

Fread and Smith (1978) proposed a methodology for estimating the pa- 
rameter “n” as a function of stage or discharge and applied the methodology 
to a hypothetical ideal single reach of a channel. This method uses modified 
Newton Raphson algorithm to improve the initial trial values, such that the 
difference between observed and computed stages is minimum. The computed 
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stages are obtained from a four point implicit finite difference solution of un- 
steady one dimensional open channel flow equations subject to upstream and 
downstream boundary conditions of observed discharge and stage hydrograph 
respectively. Fread and Smith (1978) applied the above algorithm sequen- 
tially to a multiple reach river system by specifying the computed discharge 
at the downstream end of m th reach as the upstream boundary condition for 
(m + l) th reach. The model proposed by Fread and Smith is not general in 
the sense that it can be applied to only dendritic river systems and stage ob- 
servations at m stations are absolutely necessary to estimate the “n” values 
for m reaches. Moreover, numerical errors propagate easily in any sequential 
algorithm. Therefore, any errors of estimation in the upstream reaches may 
significantly affect the estimated values in the downstream reaches. 

Recently, Wasantha Lai (1995) used Singular Value Decomposition method 
to estimate the Mannings roughness coefficient in one dimensional unsteady 
flow model. The method is used to solve for parameters after formulating the 
calibration problem as a generalized linear inverse problem. The calibration 
was repeated with different groups to determine variation of output error and 
uncertainty of parameters with parameter dimension. This method also uses 
influence coefficient approach for predicting the values of parameters at next 
iteration. But it solves underdetermined determined ,over determined and 
mixed dete- rmined problems. 

In all these previous works, the parameters it is required to solve the 
governing equations in a separate solver outside the optimization model or 
algorithm.This may necessitate running of the solver a large number of times 
iteratively in order to determine the objective function value and the influence 
coefficient matrix. In the iterative method, when the simulation is outside the 
optimization model, convergence to an optimum may prove to be difficult if 
a large number of parameters are to identified. The aim of the present study 
is to estimate the parameter “n” by directly embedding the finite difference 
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approximations of the system of governing partial differential equations 
for unsteady open channel flow into the non-linear optimization model. This 
study involves an extension of the methodology proposed by A. Narayana(199C). 

A non-linear optimization technique is used since the momentum equation and 
the objective function are non-linear. Priessmann implicit scheme (Chaudhry,1993) 
is used for discretization of the governing equations to formulate the con- 
straints. This helps in taking large computational time steps. Moreover, sta- 
bility and convergence criteria need not be checked. The main advantage of 
using an embedded optimization model for solving the inverse problem is that 
an iterative procedure is not necessary and the flow hydraulics is simulated 
within the optimization model. Also, the embedding approach results in an 
elegant formulation and has the potential of incorporating very complex flow 
systems and boundary conditions. 

1.3 Objective of the Study 

Objectives of the present study are : 

• To evaluate the proposed methodology for a number of specified channel 
systems including a single reach channel, a multiple reach channel and a 
multiple channel system; 

• To study the effect of discretizing the governing equations on the perfor- 
mance of optimizer; and 

• To evaluate performance of the proposed methodology for different sce- 
narios of data availability and noise in the flow measurement data. 
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1.4 Organization of the Thesis 


This thesis is organized as follows : 


• Chapter 1 contains the introduction to this study and discussion to rel- 
evant literature; 

• Chapter 2 presents the methodology and also discusses the simulation 
and Optimization algorithms used; 

• Solution results and performance evaluations are discussed in Chapter 3. 
and 


• Summary and Conclusions are presented in Chapter 4. 



Chapter 2 


Methodology 


2.1 Methodology 

The proposed model for estimating “n” involves the use of a nonlinear opti- 
mization technique with the finite difference form of the governing equations 
of flow forming the equality constraints. Governing equations are a set of non- 
linear hyperbolic partial differential equations representing the conservation 
of mass and momentum. These equations are discretized using the Priess- 
mann implicit finite difference scheme(Chaudhry, 1993). This scheme is used 
in the present study because large computational time steps can be consid- 
ered. The same scheme is applied to the simulation model but the momentum 
equation is linearized. The optimization algorithm used in this study is the 
Sequential Quadratic Programming (SQP) which is available in the NAG li- 
brary (NAG,1990) as the subroutine E04UCF. The routine is designed to min- 
imize an arbitary smooth nonlinear objective function subject to constraints 
which may be simple bounds on variables, linear constraints and nonlinear 
constraints. The details of the methodology are presented in this chapter. 
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2.2 Governing Equations 


In this study, one dimensional shallow water flow equations for a prismatic, 
wide rectangular channel are used as governing equations. However, it is sim- 
ple to extend the methodology to flows in irregular channels. The governing 
equations represent the mass and momentum conservation and are written as 
(Chaudhry 1993) as follows. 


Continuity equation: 


dh dq 

m + di ~ 0 

(2.1) 

Momentum equation: 


^ + + 9h2 ^ = 9h( ' S ° ~ Sf ^ 

(2.2) 


where h(x, t ) = flow depth (m), q(x, t ) =discharge per unit width of channel(m 2 /s), 
s 0 =bottom slope of channel, s/ =friction slope, x=distance along flow direc- 
tion^), and t =time(s). 

The friction slope, Sf is given by the Mannings equation: 


s f 


q 2 n 2 
flW 3 


(2.3) 


where, n is the Manning’s roughness coefficient 


2.3 Simulation Model 


A simulation model is used in the present study to simulate the observed 
flow data, in the absence of flow measurement data. The governing equations 
are linearized in the simulator. Momentum equation in nonconservation form 
along with the continuity equation(2.1) is used for the purpose. 


dq qdh q dq 
dt h dt + hdx 


q 2 dh 
h? dx 


ghs 0 + gh 


dh 

dx 


oWj = n 

lz 1 


(2.4) 
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where k is the conveyance. The derivatives are now discretized using the Priess- 
mann scheme. The variables and their time derivatives and space derivatives 
are discretized as given below. 

df (/*+’+ + 
dt 2A t 

ef ... ou£i - in t (i - ojuu - ij) 

dx Ax Ax 

f = \ + tf +1 ) + (i - «)(&, + //)) 

where, Ax = Grid spacing, At — Computational time step, 9 = weighting 
coefficient, the subscript j refers to the grid point in space and the superscript 
k refers to the time level. In Eqs.(2.-5) and (2.6), f refers to both q and h. 
In equation (2.7) f refers to variables sj,h and q. The scheme is stable if 
0.5 < 9 < 1.0 . In our study, 9 = 0.67 is used for single reach channel problem 
and 9 = 0.8 for multiple reach problems. Substituting the Eqs.(2.5)-(2.7) in 
Eqs.(2.1) and (2.4), Linearizing the equations by neglecting the second order 
terms in the power series expansion, and applying /j 1 = /* + A / 
one can obtain continuity and momentum equations in the following form: 

HjAh j+l + bjAqj+i = CjLhj + D 3 Aq 3 + G 3 (2.8) 

H' 3 Ahj + i + 6'Aq J+ i = C'jAhj + D' 3 Aq 3 + G' (2.9) 

The Eqs.(2.8) and (2.9 ) form a set of linear equations in the unknowns Ah and 
A q at all the nodes j = 1 to N (N = Maximum number of nodes). The values 
of Hj, b 3 , C 3 , D 3 , Gj, H' p b'pCj, D' and G' can be computed from the values of q 
and h at the known time level k. These values of values are known either from 
a previous computation or from specified initial conditions. Equations(2.8) 
and (2.9) along with the linearized form of boundary conditions are solved 


(2.5) 

( 2 . 6 ) 

(2.7) 
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simultaneously to obtain Ah and Aq at all the nodes j = 1 to N. Matrix 
inversion methods can be applied for solving the Equations(2.8) and (2.9). 
However, an efficient Double Sweep Algorithm as devised by Priessmann and 
Cunje is used for simulating observed data in case of single channels. In case 
of multiple channel systems, this method does not work and hence matrix 
inversion is used. This method works on the fact that if there exists a linear 
relationship of the form 

A qj = EjAhj + Fj (2.10) 

at a node j, then an analogous linear relationship 

^■Qj + i = Ej+iAhj+i + Fj+i ( 2 - 11 ) 

also exists at point j+1. By substituting the value of A qj from Eq.(2.10) into 
Eq.(2.8) we get, 

Ahj = LjAh j+l + MjAq j+l + Nj (2.12) 

where, 

Li = - ; Mi = - 6j - A r ,- 

Cj + DjEj Cj + DjEj 

Substituting the value of Aqj from Eq.(2.10) and Ahj from Eq.(2.12) 
in Eq.(2.9), gives the following relation between Ag J+1 and A hj+i can be ob- 
tained. 

+ CjE,) - gj(Cj + DjEj) 

A,)+1 (.,(0, + £>,£,) -i.,(C' + D'B,) )+ ’ 

. (Cj + + DjEj) - (G, + DjFjXq + DjEj) 

Vj(c, + DjEj) - Mq + &&) 1 

which is of the form of Eq.(2.11). All the coefficients of Eqs.(2.8) and (2.9) can 
be evaluated, as they are in terms of flow variables at known time level. In this 
study an inflow hydrograph of the form q* +l = fu(t ) is applied is applied as the 


Gj + DjE) 
Cj + DjEj 


(2.13) 
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Compute E j from B.C at 
point j=l 


Compute coefficients L:, M 

J j* 

N . as defined by Eqn(2. 13) 


Compute coefficients Ej and 
F j as defined by Eqn(2.14) 


Store L ;,M;, N z , 


]SJ n IS 

C/t j-t-1 = Nj 


Compute Ah N from Boundary 
condition at point j+1 = N 

' 


Compute A q 
k+l k 

h N =h N + Ah f 

= EA bxr + F. 

N N J 

k+l k A 

— — 



l 

Compute A bj = 

= L i Ah j + l + 



Compute Aq. = Ej ^bj + F. 
Hence calculate <j k+1 . fct 1 


k = k+l 


1 

2 * i 

, /TT\ . J 

Is 

J =J-1 

*• (NoT 

j = 1 

V. , 


Figure 2.1: Flow Chart for Double Sweep Method of Solution 
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boundary condition at the upstream point. Here ju{t) is the specified inflow 
hydrograph. Therefore, E\ becomes 0 and F\ == A q. Uniform flow condition is 
used as the boundary condition at the downstream end. The linearized form 
of this boundary condition is given as, A q = (5/3)i(/i N ) 2/3 A/iyS;. Starting 
from upstream end, the values of Ej and Fj at all the grid points upto the last 
point is obtained. From the relations at the downstream end the value of Ah 
and hence h at next time level can be obtained. The value of q at next time 
level can be obtained from the Eq.(2.10). A backward procedure is adopted 
to get the values of flow variables at all the upstream nodes at the next time 
level. The above procedure is summarized in the Fig. 2.1. 


2.4 Formulation of the Inverse Problem 

The inverse problem is formulated as a non-linear optimization model with 
roughness coefficient (n) as an unknown variable. The objective function is 
defined as the sum of squares of the differences between the estimated and 
measured values of flow depth and/or discharge. Finite difference analogs 
of the governing Eqs.(2.1) and (2.2) , and specified boundary conditions at 
the upstream and downstream ends constitute the equality constraints of the 
optimization model. The continuity equation for each node with the finite dif- 
ference node together with the upstream boundary condition at the first node 
form the linear equality constraints. The momentum equation for each node 
along with the downstream boundary condition form the non-linear equality 
constraints. Formulation of the inverse problem for a single reach channel can 
be summarised as follows. 
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2.4.1 Objective function 

Minimize : 

X [( x (°) ~ x(m))/x(m)] 2 (2.15) 

Where, x(o) are the simulated depth or discharge values at an obser- 
vation station and x(m ) are the observed depth or discharge values at the 
corresponding points. 


2.4.2 Linear constraints 


Finite difference analog of the continuity equation using the Priessmann scheme 
gives the linear constraints in the optimization model. 


(.hj+i + ht l ) ~ (/4l + ft?) . 

2Af Ax 




Ax 


0 (2.16) 


for i = 1 to N-l ; k = 1 to p 


where N = total number of grid points , p = total number of time steps 
in the computational domain and 9 = 0.67. The specified upstream boundary 
condition is also included as a linear constraint. This is given as 

q* +l = Mt) (2.17) 


where, fu(t) = Specified inflow hydrograph. 

2.4.3 Non linear constraints 


Finite difference analog of the momentum equation using the Priessmann 


scheme forms the set of nonlinear constraints. 


+ (9/ Ax) (W 


hi n 


(«? +1 n 

hi» ) 


+ ((1 - 6)/ Ax) 
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(Ail Al) 

hi 1 hi ) 


+ (s#/2Ax)((A*+*) 2 - (fc*«) 2 ) + ( 9 (1 - 9)/2Ai)((ft‘. 1 ) 2 - (hi) 2 ) 


+(g9n 2 /2) 


( (Ml , (rf +1 ) 2 'i 
U0 7/3 (h1 +, ) r ' 3 J 


(Jig? + A‘) + (1 - 0)/2(/>‘ +1 + hi)) 


+ 9 ((1 - 0)n 2 /2) 
for i = 1 to N-l ; k = lto p. 


(Ail , fa-n 

(Af +1 ) 7 / 3 + (ftf)w; 


= 0 


(2-18) 


A rating curve if available, can be specified as the boundary condition at 
the downstream end. In all the studies reported here , a uniform flow condition 
is assumed as the downstream boundary condition. This is given as 


,*+, _ (C'Pv^ 

' AT 


(2.19) 


Eq.(2.19) is also included in the set of nonlinear constraints. Formulation of 
the parameter estimation model for multiple reach system is exactly same as 
the prior described model for a single reach channel, except for the fact that 
the number of decision variables will increase as the number of “n” values to 
be estimated are multiple. Formulation of the parameter estimation model 
for a channel system is very similar to the prior described model for a single 
reach channel with the compatibility conditions at junctions forming addi- 
tional equality constraints. Mass balance and energy balance conditions at the 
junctions constitute these compatibility conditions. It should be noted here 
that other inequality constraints which incorporate lower and upper bounds 
on parameters can be also included in the optimization model. The nonlinear 
optimization problem described in this section is solved using the Sequential 
Quadratic Programming (SQP) algorithm (Powell 1974) as coded in the NAG 
library (Numerical Algorithm Group , 1990). 
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2.4.4 Optimization algorithm 


The nonlinear optimization model is solved using the sequential quadratic 
programming(SQP) algorithm, as coded in the NAG library(NAG,1990). The 
salient features of this algorithm as described in the NAG Users Manual is 
presented here. 


The SQP algorithm is used to solve the general optimization problem 
which can be stated as follows : 


Min :F(x) Subject to l < < 


Al 

C(x) 


< u 


( 2 . 20 ) 


where F(x) is the nonlinear objective function, A L is an rii x n coefficient 
matrix, C(x) is an n N element vector of nonlinear constraints. The objective 
function and constraint functions are assumed to be smooth i.e., at least twice- 
continuosly differentiable. Here upper and lower bounds are specified for all 
variables and constraints. An equality constraint can be specified by setting 
li = Ui. If certain bounds are not present, the associated elements of 1 or u can 
be set to special values that will be treated as ±oo. The initial estimates of the 
solution to (2.20), together with subroutines that define F(x), c(x) and as many 
first partial derivatives as possible must be specified; unspecified derivatives 
will be approximated by finite difference at a non-trivial expense. In this study, 
all the first partial derivatives of nonlinear constraints and objective function 
are specified by the program in the form of Jacobian and Objective gradient, 
respectively. Since the initial estimates may not be close to actual ones Cold 
Start option is used. 

At the solution of Eq.(2.20), some of the constraints will be active, i.e., 
satisfied exactly. An active simple bound constraint implies that the corre- 
sponding variable is fixed at its bound, and hence the variables are partitioned 
into fixed and free variables. A point x is a first-order Kuhn-Tucker point for 
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Eq.(2.20) if the following conditions hold(Powell, 1974): 

1. a; is feasible; 

2. there exists vectors £ and A (the Langrange multiplier vectors for bound 
and general constraints) such that 

g = C T A + f (2.21) 


where g is the gradient of F evaluated at x, and & = 0 if the j th variable is free. 

3. The Langrange multiplier corresponding to an equality constraint active 
at its lower bound must be non-negative, and non-positive for an inequality 
constraint active at its upper bound. 

There are several ways of organizing the matrix calculations of active set 
algorithm for quadratic programming. By working with the upper triangular 
matrix method for quadratic programming, it is possible to carry out matrix 
operations economically. The advantage of this approach is that there is sub- 
stantial saving of computer storage, because the orthogonal matrix is not re- 
quired. All these features are included in the E04UCF subroutine(NAG,1990), 
which is used in this study. The basic structure of E04UCF involves major 
and minor iterations. The major iterations generate sequence of x^ that con- 
verge to a;*, a first order Kuhn-Tucker Point of Eq.(2.20). At a typical major 
iterations, the new iterate x is defined by 

x = x + ap (2.22) 

where x is the current iterete, the non-negative scalar a is the step length 
and p is the search direction. Also associated with each major iteration are 
estimates of Langrange multipliers and prediction of active set. The search 
direction p is the solution of the quadratic programming subproblem of the 

Minimize: g T p + -p T Hp 


form 



2.4 Formulation of the Inverse Problem 


18 


subject to l < 


< 


V 

A L p 


> < u 


l a nP J 


(2.23) 


where g is the projected gradient of aF at x, the matrix H is the positive 
definite quasi-Newton approximation to the Hessian of the Langrangian func- 
tion and A n is the jacobian matrix of c evaluated at x. Let l in Eq.(2.20) be 
partitioned into three sections: IbJl &nd Iff, corresponding to bound, linear 
and nonlinear constraints. The vector I is similarly partitioned and defined as 
Ib = h ~ x, — A l x, and Iff = l N - c, 

where c is the vector of nonlinear constraints evaluated at x. The vector u is 
defined in an analogous fashion. 


The estimated Langrangian multipliers at each iteration are Langrange 
multipliers from subproblem (2.23). Since solving a quadratic subprogram is 
an iterative procedure, the minor iterations of the main problem form the 
iterations of the quadratic programming subproblem. 


In general, a quadratic program must be solved by iteration . Let p denote 
the current estimate of the solution of (2.23); the new iterate p is defined by 

p = p + od (2.24) 

where, as in (2.22), a is the non-negative step length and d is a search direction. 

At the beginning of each iteration of E04UCF, a working set is defined 
of constraints (general and bound) that are satisfied exactly. The vector d is 
then constructed so that the values of constraints in the working set remain 
unaltered for any move along d. For bound constraint in the working set this 
property is achieved by setting the corresponding component of d to zero, i.e. 
by fixing the variable at its bound. The subscripts ’FX’ and ’FR’ denote the 
selection of components associated with fixed and free variables. The general 
constraints in the working set will remain unaltered if 


— 0 > 


(2.25) 
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which implies 

d F R = Zd z (2.26) 

some vector d z , C is the submatrix of rows of ^ ^ ^ corresponding to general 

constraints in the working set and Z is the matrix associated with the TQ 
factorisation of C F r ■ Each change in working set leads to a simple change in 
C F r. A row or a column of C F r changes if the status of a general constraint 
or bound constraint are altered respectively. 

After obtaining the search direction from the Quadratic programming 
Subproblem, each major iteration procedes by determining a step length a 
which leads to a sufficient decrease in the Langrangian merit function given as 
under : 

L(x, A, s) = F(x) -J2 x i fe(®) ~ ^ ) + \ Y, Pi ((ciO) - Sif (2.27) 

i ^ i 

where x, A and s vary during linesearch. The summation terms in the above 
function involve only the nonlinear constraints. The vector A is an estimate 
of the Langrange multipliers for the nonlinear constraints of Eq. (2.20) and St- 
are the non-negative slack variables for nonlinear inequality constraints. The 
solution of QP subproblem (2.23) provides a vector triple that serves as the 
search direction for the three sets of the variables. 

At the end of each major iteration, a new Hessian approximation H is 
defined as a rank-two modification of H. 

«= H -Tnrs hssTH+ h vyT (2 - 28) 

where s = x-x (change in x) In E04UCF, H is required to be positive definite. 
If H is positive definite, H difined by Eq(2.28) will be positive-definite if and 
only if y F is positive. Ideally, y in Eq(2.28) would be taken as y^, the change 
in gradient of the Langrangian function 
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VL = g~ A N fj, N - AjffjL N (2.29) 

where denotes the QP multipliers associated with the nonlinear constraints 
of the original problem. If y £ is not suficiently positive, an attempt is made 
to perform the original update with vector of the form 

rriN 

V = Vl + Y, UJ i{a i xc i (x) - Oi{x)d(x )) (2.30) 

i= 1 

where LOi > 0. if no such vector can be found, the update is performed with 
scaled yx,. 

To summarise, each major iteration the routine performs before reaching 
an optimal point or otherwise includes, 

1. the solution of the quadratic programming subproblem, 

2. a line search with an augmented Langrangian merit function and 

3. a quasi-Newton update of the approximate Hessian of the Langrangian 
function. 


Chapter 3 


Results and Discussion 


3.1 Performance Evaluation 

Performance of the proposed optimization model for the estimation of rough- 
ness coefficient is evaluated using illustrative example problems for hypotheti- 
cal open channel systems. These example problems include (1) flow in a single 
channel with a single value of “n” (2) a single straight channel with multiple 
reaches and multiple “n” values and (3) a simple dendritic system of three chan- 
nels with a multiple roughness values corresponding to different reaches. The 
observation data for these cases are simulated by the solving governing Eqns. 
(2.1) and (2.2) for assumed true values of n using the Priessmann scheme. 
These simulated observation data for discharge and flow depth are then used 
in the optimization model to estimate the roughness coefficients. Identical 
initial and boundary conditions are applied while obtaining the simulated ob- 
servation data and while solving the optimization model. The advantage of 
using simulated observation data is that it provides a means to delineate the 
actual deviation between the true values and the estimated values of “n” . In 
this approach the simulated measurement data are free of measurement noises. 
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3.1.1 Single channel problem 

A wide rectangular channel with a constant bed slope of 0.0004 is considered 
in this case. The length of the channel is 40,000m. A uniform flow depth, h 
= 2.5m and corresponding uniform flow discharge for n = 0.023 are specified 
as the initial steady state conditions. A Log-Pearson Type III hydrograph as 
shown in Fig. 3.1 is specified as the upstream boundary condition. In this 
hydrograph, q b = initial discharge=3.32m 2 /s, q p = 6.0m 2 /s, t p = time to peak 
= 2hr, t g = time to centroid of hydrograph = 2.5hr and t b = time base of 
hydrograph = 6hr. Variation of q with time for this hydrograph is given by 
the following equation 

Q = Qb + (<? P - (t/t p ) <*»-*»>> (3.1) 


Flow measurement data for this case are simulated by numerically solving 
the governing equations (2.1) and (2.2) using the Priessmann scheme. These 
measurement data for this case are obtained at half hour intervals (At = 0.5 hr) 
at measurement locations thousand metres apart (As = 1000 m). These sim- 
ulated observation data for q and h without any observation noise are used in 
the optimization model to estimate the roughness coefficient. A larger grid size 
is considered in the optimization model for reducing the computational costs 
even though measurement data are available at smaller grid sizes. Manning’s 
roughness coefficient, “n” is estimated for different scenarios as presented in 
Table 3.1 in order to study the grid size effect on the parameter estimation . 
The initial estimate of “n” which is required to start the optimization model 
is specified to be 0.01 in all the runs. It can be readily seen from Table 3.1 
that the optimization model performs very well for the case of a single channel 
even if Ax and At are large. The maximum difference between the estimated 
value and the true value is only 0.0001 . An average of 10 iterations of the Op- 
timization model are required to arrive at these optimal n values. Estimated 
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Figure 3.1: Log Pearson Typelll Hydrograph 


and observed stage hydrographs at x = 20 km for Scenario 1 are compared 
in Fig. 3.2. This comparision indicates satisfactory performance of the opti- 
mization model. The Scenario 1 has 253 variables, 126 linear and non- linear 
constraints each, respectively. So, the CPU time for this problem is 3min and 
12sec. When the grid spacing is increased, as in other Scenarios of Table 3.1 
the CPU time reduces and it comes down to a few seconds for Scenario 3. 

The effect of different initial estimates of “n” on the obtained optimum 
solution for “n” is shown in Table 3.2. Results shown in Table 3.2 are obtained 
using Ax = 2000 m and At = 1 hr in the optimization model. It can be seen 
from Table 3.2 that the optimization model performs satisfactorily even when 
the initial estimate is far off from the true value. However, as may be expected, 
less number of iterations are required when the initial estimate is close to the 
true value. In all the cases, since the order magnitude of flow variables are 
known at all points from simulated observation data, the intial estimates for 
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Figure 3.2: Estimated and observed stage hydrographs at X= 20 km (Scenario 


Table 3.1: Grid size effect on the parameter estimation 


Scenario no. 

Grid size in the Optimizer 

Estimated Value 
of V 

No. of 
iterations 

Objective 

function 

Ax (m) 

At (hr) 

1 

2000 

1.0 

0.0229 

12 

0.236 

2 

4000 

1.0 

0.0229 

9 

0.118 

3 

8000 

1.0 

0.0229 

9 

0.053 

4 

4000 

0.5 

0.0228 

9 

0.04 

5 

8000 

0.5 

0.0228 

9 

0.044 


flow variables are taken as equal to the observed values. The solution results 
are identical even the initial estimates for the flow variables are 80% of observed 
values . 
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Table 3.2: Effect of initial estimate on the performance of the model 


Scenario no. 

Initial estimate 
of V 

Optimum 
value of “n” 

No. of 
Iterations 

6 

0.2 

6.0233 

42 

7 

0.015 

0.0229 

7 

8 

0.005 

0.0229 

15 


3.1.2 Multiple reach problem 

3. 1.2.1 single channel 


A straight channel consisting of six reaches of equal length, each 6000m and 
having different roughness coefficients is considered in this case. All the reaches 
have identical bottom slopes, s 0 = 0.0004. The flow measurement data for 
this system is simulated using the inflow hydrograph shown in Fig. 3.1 as the 
upstream boundary condition. The true values of “n” used in these simulations 
are 0.018, 0.023 0.028, 0.025, 0.031 and 0.035, respectively. The observation 
data is simulated using Ax = 1000 m, At = 0.5hr and 6 = 0.67 in the 
simulator. The initial steady state conditions are obtained by solving the 
gradually varied flow equations corresponding to a discharge of 3.53m 2 /s and a 
downstream depth of 2.98m. The simulated flow measurement data are used in 
the objective function of the optimization model to estimate the “n” values for 
all the six roaches. The initial estimates for all the six “n” values are specified 
to be 0.022. Table 3.3 shows the effect of grid size on performance of the 
optimization model. It can be observed from Table 3.3 that the optimization 
model performs satisfactorily even in the case of multiple readi problem. As 
may be expected, errors in the estimated values of n increase as At and 
increased. For example, the maximum error in the estimation is only 0.0017 
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for the scenario 9. Estimated and observed stage hydrograph at x = 18 km 
are compared in the Fig. 3.3. The CPU time required for Scenario No. 9 is 
around 4 min and decreases as the number of variables decease. An average 
of 15 iterations are required to arrive at these optimal values. 

Table 3.3: Grid size effect:Multiple reach single channel problem 


Scenario no. 

Grid Size 

Estimated Value of “n” 

Ax (m) 

At (hr) 

Reach no. 

1 

2 

3 

4 

5 

6 

9 

2000 

1.0 

0.0178 

0.0220 

0.0263 

0.0244 

0.0298 

0.0352 

10 

■illM 


0.0176 

0.0212 

0.0245 

0.0232 

BliESEl 

liliEEEl 

11 

iiaiM 

2.0 

0.0176 

0.0221 

0.0274 

0.0254 

IflBSCT 

BliEEEl 



Time in hours 

Figure 3.3: Estimated and observed stage hydrographs at the end of third 
reach (Scenario 9) 
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3.1. 2. 2 Multiple channel system 

A dendritic system of three wide rectangular channels as shown in Fig. 3.4 
is considered in this case. Each channel is divided into two reaches of equal 
length, with each reach having a different roughness coefficient. All three 
channels are of the same length = 12,000m, with identical bed slope, s 0 — 
0.0004. The flow measurement data for this system is simulated also by using 
the inflow hydrograph shown in Fig. 3.1 as the upstream boundary condition. 
The true values of “n” used in these simulations are 0.018, 0.023, 0.028, 0.025, 
0.031 and 0.035 for the reaches 1,2, 3, 4, 5 and 6, respectively. The observation 
data is simulated using Ax = 1000 m and A f = 0.5 hr and 6 = 0.8 in the 
simulator. The initial steady state conditions for this system are obtained 
by solving the gradually varied flow' equation corresponding to discharges of 
3.53 m 7 /s , 2.27 m 2 /s and 5.80m 2 /s in channels 1,2 and 3 respectively. The 
initial downstream depth in channel 3 in these computations is 4.02 m. 

The simulated flow' measurement data are used in the objective function 
of the optimization model to estimate the “n” values for all the six reaches. 
The initial estimates for all the six “n” values are specified to be 0.022. Table 
3.4 shows the effect of the grid size on the performance of the optimization 
model. It can be observed from Table 3.4 that the optimization model per- 
forms satisfactorily even in the case of multiple channel problem. As may be 
expected, errors in the estimated values of n increase as A t and Ax are in- 
creased. For example, the maximum error in the estimation is only 0.0013 for 
the Scenario 12 (Ax = 2000 m and At = 1 hr) as compared to an error of 
0.0025 for scenario 13 (Ax = 6000 m , At = 1 hr). Identical results are ob- 
tained when the initial estimates for ll n” are taken equal to 0.01 indicating the 
robustness of the proposed model. Estimated and observed stage hydrographs 
6000m dowmstream of the junction point, along the channel 3 for Scenario 
12 are compared in Fig. 3.5. Satisfactory performance of the optimization 
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model can be observed again from this figure. Maximum difference between 
the estimated and the observed flow depths is only 0.11 m. 


Table 3.4: Grid size effect: Multiple channel problem 
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Flow Direction 


Reach number 


Channel number 


Figure 3.4: Schematic layout of the multiple channel system 



Time in hours 

Figure 3.5: Estimated and observed stage hydrographs 6000m downstream of 
junction point in channel 3 (Scenario 12) 
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3. 1.2. 3 Performance evaluation of the model for missing data 

In many cases, observation data may be available either for flow depths or 
discharges. Also, data may be available only at few selected stations and at 
few selected times. However, it may be required to use smaller Ax and At 
in the optimization model while estimating the parameters in order to reduce 
the discretization errors. Performance of the optimization model under such 
conditions is studied by considering several scenarios as presented in Table 
3.5. Ax and At in all these studies are 2000m and lhr, respectively. The 
initial estimates for the flow variables are taken as spatial mean of the avail- 
able data at that time level. Initial estimates for the “n” values are taken as 
0.022. Results presented in table 3.5 indicate that the proposed model per- 
formed satisfactorily even when the availability of the measurement data is 
sparse (Scenarios 15-18). However, there is a significant deterioration in the 
performance of the model when the flow depth measurements are not available 
(Scenario 19). The maximum error in estimation of “n” for this scenario being 
0.010, in reach 5. Results are also not satisfactory, as expected, in the case 
of an underdetermined system (Scenario 20) i.e., the number of observation 
stations is less than the number of parameters to be estimated. The maximum 
error in estimation of “n” being 0.012 for this scenario in reach 5. 
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Table 3.5: Performance evaluation of the model for the multiple channel prob- 
lem 


Scenario no. 

Availability of Observed Data 

Estimated “n” Values 

15 

(i) both q and h 

(ii) t = 1,3 & 5 hrs 

(iii) All Spatial Grid Points 

0.0178, 0.0206, 
0.0296, 0.0246, 
0.0299 and 0.0359 

16 

(i) both q and h 

(ii) t = 1 & 4 hrs 

(iii) All spatial Grid Points 

0.0182, 0.021, 
0.0293, 0.0249, 
0.0307 and 0.0361 

17 

(i) both q and h 

(ii) all times 

(iii) only the two end points of a reach 

0.018, 0.0206, 
0.0291, 0.0217 
0.0309 and 0.0365 

18 

(i) only h 

(ii) all times 

(iii) all spatial grid points 

0.0183, 0.0236, 
0.0278, 0.0256, 
0.0336, 0.0384 

19 

(i) only q 

(ii) all times 

(iii) all spatial grid points 

0.0203, 0.0191, 
0.0367, 0.0238, 
0.0211 and 0.0363 

20 

(i) both q and h 

(ii) all times 

(iii) only for longitudinal mid points 
of reaches 1,3 & 5 

(under determined system) 

0.0191, 0.0235, 
0.0276,0.0197, 
0.0323 and 0.0429 
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3. 1.2.4 Performance evaluation of the model in the presence of ob- 
servation noise 

In many cases , flow measurement data from field conditions may contain 
observation noise. Performance of the proposed parameter estimation model 
in the presence of observation noise is evaluated here. In the illustrative case 
considered here , all the input data except the flow measurement data are 
identical to those for Scenario 12. Simulated noise free measurement data for 
Scenario 12 is modified by adding a random error term as follows 

x'i = Xi + €i (3.2) 

where , X{ =simulated noise free measurement of ith data point, and e,- = a 
random error term sampled from a uniform distribution of zero mean , and 
upper and lower limits of -i -aXi and —axi , respectively. The NAG routine 
G05FAF is used to generate the random numbers. This program generates a 
vector of pseudo random numbers uniformly distributed over an interval (a,b). 
In this evaluation , five sets of noisy flow measurement data are generated using 
the prior procedure. These data sets are then incorporated into the objective 
function of the parameter estimation model to estimate the roughness values. 
The “n” values obtained for every random number in cases of induced errors 
of 10% and 20% are given in Tables 3.6 and 3.7, respectively. The average 
of “n” values obtained as solutions for the five sets of noisy data with a = 
0.1 and 0.2 are also presented. The estimated average values of n for both 
the cases of a = 0.1 and 0.2 are close to those for a = 0.0 representing, error 
free measurement. The sum of squares of the estimation errors is presented 
as a function of a in Fig. 3.6. Tables(3.6 and 3.7) and Fig. 3.6 indicate the 
potential applicability of the proposed model for real life situations where the 
measurement data are erroneous. 
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Table 3.6: Estimated values of “n” for an a. = 0.1 


Erroneous 
measurement 
data set no. 

Estimated Value of “n” 

Reach no. 

1 

2 

3 

4 

5 

6 

1 

0.0177 

0.0216 

0.0307 

0.0249 

0.0296 

0.0359 

2 

0.0174 

0.0211 

0.0286 

0.0247 


0.0362 

3 

0.0170 

0.0215 



0.0301 

BfiBSEl 

4 



0.0308 

0.0259 

0.0302 


5 

0.0194 

0.0218 

0.0291 



■tltKfflil 

Mean values 




0.0252 

0.0306 

0.0362 

True values 


0.023 

0.028 

0.025 

0.031 

0.035 


Table 3.7: Estimated values of “n” for a = 0.2 


Erroneous 
measurement 
data set no. 

Estimated Value of “n” 

Reach no. 

1 

2 

3 

4 

5 

6 

1 

liiiilwfel 

0.0216 

0.0332 


iwyial 

lytigal 

2 

EH 

0.0204 

0.0285 

gft'sfelTI 


BlifliSl 

3 


KliEilai 





4 


Hi 

0.0328 


OB 


5 


liEEHg 

0.0291 

0.0275 

EB 


Mean values 

0.0177 

EEH 


0.0254 

0.0302 

OH 

True values 

0.018 


0.028 

0.025 

0.031 

0.035 































































Chapter 4 


Summary And Conclusion 


In this study, an embedded optimization model is presented for estimating 
the Manning’s roughness coefficients from unsteady flow measurement data in 
open-channel systems. Finite-difference approximations of the governing flow 
equations are embedded as equality constraints in a nonlinear optimization 
model. Priessmann finite difference method is used for discretizing the govern- 
ing partial differential equations. The nonlinear objective function based on 
the sum of squares of error is minimized by using the Sequential Quadratic Pro- 
gramming (SQP) algorithm as coded by Numerical Algorithm Group (NAG). 
Performance of the proposed model is evaluated for different scenarios of data 
availability and observation noise in the flow measurement data. It is found 
that the model performs satisfactorily in all the cases except (1) when only 
discharge measurements are available and (2) when the number of observation 
stations is less than the number of parameters to be estimated. Effects of 
grid size and initial estimates for parameters on the model performance are 
also studied. Results of these evaluations show that the performance of the 
model is satisfactory even though global optimality of the solution cannot be 
guaranteed. Also, the applicability of the proposed approach to very large 
domain problems should be further explored, considering some of the inherent 
limitations of the embedding technique. 
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